home *** CD-ROM | disk | FTP | other *** search
/ Programmer Power Tools / Programmer Power Tools.iso / c / spline.c < prev    next >
C/C++ Source or Header  |  1986-05-18  |  16KB  |  547 lines
  1. To: info-ibmpc@usc-isib.ARPA
  2. Subject: contribution to library
  3. Date: Sat, 17 May 86 17:18:28 -0500
  4. From: James R. Van Zandt <jrv@mitre-bedford.ARPA>
  5.  
  6. I'd like to submit this program to the library.
  7. -------------------------------------------
  8.  
  9. SPLINE.C    Generates splines under tension and allows general curves
  10.             and multiple independent curves in the same file.  Text
  11.             input and output files like the UNIX program.  Written in
  12.             C by J. R. Van Zandt, based on algorithms by A.  K.  Cline.
  13.  
  14. ------------------------------------------------------------
  15. NAME
  16.     spline - interpolate using splines under tension
  17.  
  18. SYNOPSIS
  19.     spline [file] [options]
  20.  
  21. DESCRIPTION
  22.     Without options, SPLINE reads pairs of numbers (x- and y-values) from
  23.     the standard input (or the given file), generates a smooth curve
  24.     through the points, and writes to the standard output points from the
  25.     smooth curve.  The curve is a spline under tension (see references),
  26.     which is somewhat "tighter" than a cubic spline, and less likely to
  27.     have spurious inflection points.  As with GRAPH, each pair of points
  28.     may optionally be followed by a comment.  If the comment is surrounded
  29.     by quotes "...", it may contain spaces.  The given points, and
  30.     their comments if any, will be included in the output.
  31.  
  32.     Input lines starting with ";" are written to the beginning of
  33.     the output file but are otherwise ignored.  Blank lines are
  34.     ignored.
  35.  
  36.     If the -c switch is not used, the input points must be from a function
  37.     - that is, the x values must be strictly increasing.  The
  38.     output points will also be from a function.  (If the -b switch
  39.     is used, this restriction applies only within each segment.)
  40.  
  41.     If the -c switch is used (indicating a general curve), the
  42.     input points need not be from a function, but each pair of
  43.     points must be separated from the previous pair by a finite
  44.     distance.  (If the -b switch is used, this restriction applies
  45.     only within each segment.)
  46.     
  47. options are:
  48.   -a  [step [start]] Input data contains only y values - generate automatic
  49.           abscissas at intervals of step (default 1) starting at start
  50.           (default 0).
  51.  
  52.   -b         break the interpolation at each label.  That is, the input
  53.         curve is divided into sections with the last point in
  54.         each section marked by a label (which may be empty:
  55.         "").  A separate interpolating curve is to be found for
  56.         each section.  In this case, the requirements on the
  57.         number of intervals (specified by the -n switch or
  58.         defaulted) and the interpolation range (specified by the 
  59.         -x switch) are applied to each section independently.
  60.  
  61.   -c            general curve rather than function.  In this case, the curve
  62.           is parameterized on the polygonal arclength from the
  63.           first to the last given point, with the whole length
  64.           scaled to be 1.  Thus, the values min and max for the
  65.           -x switch should satisfy 0 <= min < max <= 1.
  66.  
  67.   -n  num       interpolate over num intervals (default is 100)
  68.  
  69.   -t num    Specify tension in interpolating curve.  Tension of 50 gives
  70.         almost polygonal line, tension of .01 gives almost cubic
  71.         spline.  Default is 1.
  72.  
  73.   -x [min [max]]   interpolate from min to max only.  min and max should be in
  74.         the range of the given x values, except that if the -c switch
  75.         is used they should satisfy 0 <= min < max <= 1.
  76.  
  77.   -xl        take log of x values before interpolating, take exponential
  78.            afterwards (probably necessary if -xl switch is needed for
  79.            GRAPH)
  80.  
  81.   -yl        take log of y values before interpolating, take exponential
  82.            afterwards (probably necessary if -yl switch is needed for
  83.            GRAPH)
  84.  
  85. NOTES
  86.     Similar to the Unix routine, except using splines under tension,
  87.     passing labels through, and allowing general curves.
  88.  
  89. REFERENCES
  90.  
  91.     A.  K.  Cline, "Scalar- and Planar- Valued Curve Fitting Using
  92.     Splines Under Tension", Communications of the ACM v 17 n 4 p
  93.     218-223 (Apr 74).
  94.  
  95.     Schweikert, D.  G.  "An interpolation curve using a spline in
  96.     tension", J.  Math.  and Physics v 45 p 312-317 (1966).
  97.  
  98. AUTHOR
  99.     Copyright (c) 1985  James R. Van Zandt
  100.  
  101.     Resale forbidden, copying for personal use encouraged.
  102.  
  103. ------------------------------------------------------------
  104. /*    spline - interpolate points in a tabulated function or curve
  105.  
  106.     Usage...
  107.         spline [file][options]
  108.         options are:
  109.           -a  [step [start]] automatic abscissas 
  110.           -b                 break interpolation at each label
  111.           -c                 general curve 
  112.           -n  num            interpolate over num intervals 
  113.           -t  num            tension in interpolating curve 
  114.           -x  [min [max]]    interpolate from min to max 
  115.  
  116.     Notes...
  117.         This program fits a smooth curve through a given set of points,
  118.         using the splines under tension introduced by Schweikert [J. 
  119.         Math.  and Physics v 45 pp 312-317 (1966)].  They are similar
  120.         to cubic splines, but are less likely to introduce spurious
  121.         inflection points.
  122.  
  123.     History...
  124.         21 Sep 85    Added -xl and -yl switches for taking logs
  125.         23 Nov 85    Passing lines starting with ';' unchanged, otherwise
  126.                     ignoring them.
  127.  
  128.     Author...
  129.         Copyright (c) 1985  James R. Van Zandt
  130.  
  131.         All rights reserved.
  132.         This program may be copied for personal, non-profit use only.
  133.  
  134.         Based on algorithms by A.  K.  Cline [Communications of the ACM
  135.         v 17 n 4 pp 218-223 (Apr 74)].
  136.  
  137. */
  138.  
  139. #include <stdio.h>
  140. #include <math.h>
  141.  
  142. #define VERSION "1.2"
  143.  
  144. #define ENTRIES 200
  145. #define MAXLABELS 50
  146.  
  147. double *x, *y, *temp, *xp, *yp, *path;
  148. double slp1=0.,slpn=0.,sigma=-1.;
  149. double abscissa=0.;
  150. double abscissa_step=1.;
  151. double xmin=0.;
  152. double xmax=0.;
  153. double curv2(), mylog();
  154.  
  155. int xlog=0;        /* nonzero if taking log of x values */
  156. int ylog=0;        /* nonzero if taking log of y values */
  157. int debugging=0;
  158. int breaking=0;        /* nonzero if breaking at labels */
  159. int labels=0;        /* number of labels in input data */
  160. int automatic_abscissas=0;
  161. int abscissa_arguments=0;
  162. int curve=0;        /* nonzero if general curve permitted */
  163. int x_arguments=0;
  164. int n;                /* number of entries in x, y, yp, and temp */
  165. int index_array[MAXLABELS];    /* indexes into x and y */
  166. int *p_data=index_array;
  167. int total=100;
  168.  
  169. FILE file;
  170.  
  171. char *label_array[MAXLABELS];
  172. char **p_text=label_array;
  173.  
  174. main(argc,argv) int argc; char **argv;
  175. {    int nn,origin;
  176.     read_data(argc,argv);
  177.     if(breaking && labels)
  178.         {origin=0;
  179.         while(labels--)
  180.             {p_data[0] -= origin;
  181.             nn=p_data[0]+1;
  182.             if(nn) curv0(nn,total);
  183.             origin += nn;
  184.             path += nn;
  185.             x += nn;
  186.             y += nn;
  187.             p_data++;
  188.             p_text++;
  189.             }
  190.         }
  191.     else
  192.         curv0(n,total);
  193. }
  194.  
  195. curv0(n,requested)
  196. int n,requested;
  197. {    int i, j, each, stop, seg=0, regressing=0;
  198.     double s,ds,xx,yy;
  199.  
  200.     for(i=1; i<n; i++) if(path[i]<=path[i-1]) regressing++;
  201.     if (regressing || (curve && (xmin<0. || xmax>1.)))
  202.         {              /* error: echo input */
  203.         for (i=0; i<n; i++)
  204.             {xx=x[i]; if(xlog) xx=exp(xx);
  205.             yy=y[i]; if(ylog) yy=exp(yy);
  206.             printf("%g %g ",xx,yy);
  207.             if(i == p_data[seg]) puts(p_text[seg++]);
  208.             putchar('\n');
  209.             }
  210.         return;
  211.         }
  212.     if(curve) {curv1(path,x,xp,n); curv1(path,y,yp,n);}
  213.     else curv1(x,y,yp,n);
  214.     s=path[0];
  215.     seg=0;
  216.     if(requested<n) requested=n;
  217.     if(x_arguments==2)        /* specific range was requested */
  218.         {s=xmin;
  219.         ds=(xmax-xmin)/requested;
  220.         if(curve) {ds*=path[n-1]; s*=path[n-1];}
  221.         for(i=requested+1; i; i--)
  222.             {value(s,&xx,&yy,n);
  223.             printf("%g %g ",xx,yy);
  224.             if(j==p_data[seg]) puts(p_text[seg++]);
  225.             putchar('\n');
  226.             s+=ds;
  227.             }
  228.         }
  229.     else            /* spline for entire curve was requested */
  230.         {for (i=j=0; j<n-1; j++)
  231.             {stop=requested*(j+1)/(n-1);
  232.             ds=(path[j+1]-path[j])/(stop-i);
  233.             for ( ; i<stop; i++)
  234.                 {value(s,&xx,&yy,n);
  235.                 printf("%g %g ",xx,yy);
  236.                 if(j==p_data[seg]) puts(p_text[seg++]);
  237.                 putchar('\n');
  238.                 s+=ds;
  239.                 }
  240.             }
  241.             xx=x[n-1]; if(xlog) xx=exp(xx);
  242.             yy=y[n-1]; if(ylog) yy=exp(yy);
  243.             printf("%g %g ",xx,yy);
  244.             if(j==p_data[seg]) puts(p_text[seg++]);
  245.             putchar('\n');
  246.         }
  247. }
  248.  
  249. value(s,q,r,n) double s,*q,*r; int n;
  250. {    if(curve)
  251.         {*q=curv2(path,x,xp,s,n);
  252.         *r=curv2(path,y,yp,s,n);
  253.         }
  254.     else
  255.         {*q=s;
  256.         *r=curv2(x,y,yp,s,n);
  257.         }
  258.     if(xlog) *q=exp(*q);
  259.     if(ylog) *r=exp(*r);
  260. }
  261.  
  262. read_data(argc,argv) int argc; char **argv;
  263. {    int i,j,length;
  264.     double xx,yy,d,*pd,sum;
  265.     char *s,*t;
  266. #define BUFSIZE 200
  267.     static char buf[BUFSIZE];
  268.  
  269.     x=path=malloc(ENTRIES*sizeof(double));
  270.     y=malloc(ENTRIES*sizeof(double));
  271.     if(x==0 || y==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  272.     if(argc>1 && streq(argv[1],"?")) help();
  273.     if(argc<=1 || *argv[1]=='-') file=stdin;
  274.     else
  275.         {if(argc>1)
  276.             {file=fopen(argv[1],"r");
  277.             if(file==0) {printf("file %s not found\n",argv[1]); exit();}
  278.             argc--; argv++;
  279.             }
  280.         else help();
  281.         }
  282.     argc--; argv++;
  283.     while(argc>0)
  284.         {i=get_parameter(argc,argv);
  285.         argv=argv+i; argc=argc-i;
  286.         }
  287.     if(xlog && !curve)
  288.         {if(x_arguments>1) xmax=mylog(xmax);
  289.         if(x_arguments>=1) xmin=mylog(xmin);
  290.         }
  291.     if(automatic_abscissas && abscissa_arguments<2 && x_arguments>=1)
  292.         abscissa=xmin;
  293.     p_data[0]=-1;
  294.     i=0;
  295.     while(i<ENTRIES)
  296.         {if(fgets(buf,BUFSIZE,file)==0) break;
  297.         t=buf;
  298.         while(*t && isspace(*t)) t++;
  299.         if(*t == 0) continue;        /* skip blank lines */
  300.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  301.         if(buf[0]==';') {printf("%s\n",buf); continue;}  /* skip comment */
  302.         if(automatic_abscissas)
  303.             {x[i]=abscissa;
  304.             abscissa+=abscissa_step;
  305.             sscanf(buf,"%F",&y[i]);
  306.             if(ylog) y[i]=mylog(y[i]);
  307.             }
  308.         else
  309.             {sscanf(buf,"%F %F",&x[i],&y[i]);
  310.             if(xlog) x[i]=mylog(x[i]);
  311.             if(ylog) y[i]=mylog(y[i]);
  312.             }
  313.         s=buf;                      /* start looking for label */
  314.         while(*s==' ')s++;        /* skip first number */
  315.         while(*s && (*s!=' '))s++;
  316.         if(!automatic_abscissas)    /* skip second number */
  317.             {while(*s==' ')s++;
  318.             while(*s && (*s!=' '))s++;
  319.             }
  320.         while(*s==' ')s++;
  321.         if((length=strlen(s))&&(labels<MAXLABELS))
  322.             {if(*s=='\"')           /* label in quotes */
  323.                 {t=s+1;
  324.                 while(*t && (*t!='\"')) t++;
  325.                 t++;
  326.                 }
  327.             else                    /* label without quotes */
  328.                 {t=s;
  329.                 while(*t && (*t!=' '))t++;
  330.                 }
  331.             *t=0;
  332.             length=t-s;
  333.             p_data[labels]=i;
  334.             p_text[labels]=malloc(length+1);
  335.             if(p_text[labels]) strcpy(p_text[labels++],s);
  336.             }
  337.         i++;
  338.         }
  339.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  340.         {p_data[labels]=i-1;
  341.         if(p_text[labels]=malloc(1)) *p_text[labels++]=0;
  342.         }
  343.     n=i;
  344.     yp=malloc(n*sizeof(double));
  345.     temp=malloc(n*sizeof(double));
  346.     if(temp==0 || yp==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  347.     if(curve)
  348.         {xp=malloc(n*sizeof(double));
  349.         path=malloc(n*sizeof(double));
  350.         if(xp==0|| path==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  351.         path[0]=sum=0.;
  352.         for (i=1; i<n; i++)
  353.             {xx=x[i]-x[i-1];
  354.             yy=y[i]-y[i-1];
  355.             path[i]=(sum+=sqrt(xx*xx + yy*yy));
  356.             }
  357. /*        for(i=0; i<n; i++)
  358.             printf("path[%d]=%g  x[%d]=%g \n",i,path[i],i,x[i]); */
  359.         }
  360. }
  361.  
  362.  
  363. /* get_parameter - process one command line option
  364.         (return # parameters used) */
  365. get_parameter(argc,argv) int argc; char **argv;
  366. {    int i;
  367.     if(streq(*argv,"-a"))
  368.         {i=get_double(argc,argv,2,&abscissa_step,&abscissa,&abscissa);
  369.         abscissa_arguments=i-1;
  370.         automatic_abscissas=1;
  371.         return i;
  372.         }
  373.     else if(streq(*argv,"-b")) {breaking=1; return 1;}
  374.     else if(streq(*argv,"-c")) {curve=1; return 1;}
  375.     else if(streq(*argv,"-d")) {debugging=1; return 1;}
  376.     else if(streq(*argv,"-n"))
  377.         {if((argc>1) && numeric(argv[1])) total=atoi(argv[1]);
  378.         if(debugging)printf("%d output pairs",total);
  379.         return 2;
  380.         }
  381.     if(streq(*argv,"-t"))
  382.         {return(get_double(argc,argv,1,&sigma,&abscissa,&abscissa));
  383.         }
  384.     else if(streq(*argv,"-x"))
  385.         {i=get_double(argc,argv,2,&xmin,&xmax,&xmax);
  386.         x_arguments=i-1;
  387.         return i;
  388.         }
  389.     else if(streq(*argv,"-xl")) {xlog++; return 1;}
  390.     else if(streq(*argv,"-yl")) {ylog++; return 1;}
  391.     else gripe(argv);
  392. }
  393.  
  394. get_double(argc,argv,permitted,a,b,c)
  395. int argc,permitted; char **argv; double *a,*b,*c;
  396. {    int i=1;
  397.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  398.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  399.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  400.     return i;
  401. }
  402.  
  403. int streq(a,b) char *a,*b;
  404. {    while(*a)
  405.         {if(*a!=*b)return 0;
  406.         a++; b++;
  407.         }
  408.     return 1;
  409. }
  410.  
  411. gripe_arg(s) char *s;
  412. {    fprintf(stderr,"argument missing for switch %s",s);
  413.     help();
  414. }
  415.  
  416. gripe(argv) char **argv;
  417. {    fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
  418.     help();
  419. }
  420.  
  421. help()
  422. {    fprintf(stderr,"spline   version %s",VERSION);
  423.     fprintf(stderr,"\nusage: spline [file][options]\n");
  424.     fprintf(stderr,"options are:\n");
  425.     fprintf(stderr,"  -a  [step [start]] automatic abscissas \n");
  426.     fprintf(stderr,"  -b                 break interpolation after each label \n");
  427.     fprintf(stderr,"  -c                 general curve \n");
  428.     fprintf(stderr,"  -n  num            interpolate over num intervals \n");
  429.     fprintf(stderr,"  -t  num            tension in interpolating curve \n");
  430.     fprintf(stderr,"  -x  [min [max]]    interpolate from min to max \n");
  431.     fprintf(stderr,"  -xl                take logs of x values before interpolating \n");
  432.     fprintf(stderr,"  -yl                take logs of y values before interpolating \n");
  433.     exit();
  434. }
  435.  
  436. numeric(s) char *s;
  437. {    char c;
  438.     while(c=*s++)
  439.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  440.         return 0;
  441.         }
  442.     return 1;
  443. }
  444.  
  445. curv1(x,y,yp,n) double *x,*y,*yp; int n;
  446. {    int i;
  447.     double c1,c2,c3,deln,delnm1,delnn,dels,delx1,delx2,delx12;
  448.     double diag1,diag2,diagin,dx1,dx2,exps;
  449.     double sigmap,sinhs,sinhin,slpp1,slppn,spdiag;
  450.  
  451.     delx1=x[1] - x[0];
  452.     dx1=(y[1] - y[0])/delx1;
  453.     if(sigma >= 0.) {slpp1=slp1; slppn=slpn;}
  454.     else
  455.         {if(n!=2)
  456.             {delx2= x[2] - x[1];
  457.             delx12= x[2] - x[0];
  458.             c1= -(delx12 + delx1)/delx12/delx1;
  459.             c2= delx12/delx1/delx2;
  460.             c3= -delx1/delx12/delx2;
  461.             slpp1= c1*y[0] + c2*y[1] + c3*y[2];
  462.             deln= x[n-1] - x[n-2];
  463.             delnm1= x[n-2] - x[n-3];
  464.             delnn= x[n-1] - x[n-3];
  465.             c1= (delnn + deln)/delnn/deln;
  466.             c2= -delnn/deln/delnm1;
  467.             c3= deln/delnn/delnm1;
  468.             slppn= c3*y[n-3] + c2*y[n-2] + c1*y[n-1];
  469.             }
  470.         else yp[0]=yp[1]=0.;
  471.         }
  472.             /* denormalize tension factor */
  473.     sigmap=fabs(sigma)*(n-1)/(x[n-1]-x[0]);
  474.             /* set up right hand side and tridiagonal system for
  475.                yp and perform forward elimination    */
  476.     dels=sigmap*delx1;
  477.     exps=exp(dels);
  478.     sinhs=.5*(exps-1./exps);
  479.     sinhin=1./(delx1*sinhs);
  480.     diag1=sinhin*(dels*.5*(exps+1./exps)-sinhs);
  481.     diagin=1./diag1;
  482.     yp[0]=diagin*(dx1-slpp1);
  483.     spdiag=sinhin*(sinhs-dels);
  484.     temp[0]=diagin*spdiag;
  485.     if(n!=2)
  486.         {for(i=1; i<=n-2; i++)
  487.             {delx2=x[i+1]-x[i];
  488.             dx2=(y[i+1]-y[i])/delx2;
  489.             dels=sigmap*delx2;
  490.             exps=exp(dels);
  491.             sinhs=.5*(exps-1./exps);
  492.             sinhin=1./(delx2*sinhs);
  493.             diag2=sinhin*(dels*(.5*(exps+1./exps))-sinhs);
  494.             diagin=1./(diag1+diag2-spdiag*temp[i-1]);
  495.             yp[i]=diagin*(dx2-dx1-spdiag*yp[i-1]);
  496.             spdiag=sinhin*(sinhs-dels);
  497.             temp[i]=diagin*spdiag;
  498.             dx1=dx2;
  499.             diag1=diag2;
  500.             }
  501.         }
  502.     diagin=1./(diag1-spdiag*temp[n-2]);
  503.     yp[n-1]=diagin*(slppn-dx2-spdiag*yp[n-2]);
  504.                     /* perform back substitution */
  505.     for (i=n-2; i>=0; i--) yp[i] -= temp[i]*yp[i+1];
  506. }
  507.  
  508.  
  509. double curv2(x,y,yp,t,n) double *x,*y,*yp,t; int n;
  510. {    int i,j;
  511.     static int i1=1;
  512.     double del1,del2,dels,exps,exps1,s,sigmap,sinhd1,sinhd2,sinhs;
  513.  
  514.     s=x[n-1]-x[0];
  515.     sigmap=fabs(sigma)*(n-1)/s;
  516. #ifdef WORK
  517.     for (j=2; j; j--)    /* want: x[i-1] <= t < x[i], 0 < i <= n */
  518.         {for (i=i1; i<n; i++)
  519.             {if(x[i]>t) break;
  520.             }
  521.         if(i==n) i=n-1;
  522.         if(x[i-1]<=t || t<=x[0]) break;
  523.         i1=1;
  524.         }
  525. #endif
  526.     i=i1;
  527.     while(i<n && t>=x[i]) i++;
  528.     while(i>1 && x[i-1]>t) i--;
  529.     i1=i;
  530.     del1=t-x[i-1];
  531.     del2=x[i]-t;
  532.     dels=x[i] - x[i-1];
  533.     exps1=exp(sigmap*del1); sinhd1=.5*(exps1-1./exps1);
  534.     exps= exp(sigmap*del2); sinhd2=.5*(exps-1./exps);
  535.     exps=exps1*exps;        sinhs=.5*(exps-1./exps);
  536.     return ((yp[i]*sinhd1 + yp[i-1]*sinhd2)/sinhs +
  537.             ((y[i] - yp[i])*del1 + (y[i-1] - yp[i-1])*del2)/dels);
  538. }
  539.  
  540. double mylog(x) double x;
  541. {    if(x>0.) return log(x);
  542.     fprintf(stderr,"can%'t take log of nonpositive number");
  543.     exit(1);
  544.     return 0.;
  545. }
  546.  
  547.